#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include "dept150_gaus_mgpu.h"
#include "mkl_lapack.h"

int main(int argc, char** argv)
{
	double* A = NULL;
	int* ipiv = NULL;
	int i, j;
	int N = 6;
	int block_size = 3;
	if(argc >=  2)
	{
		N = atoi(argv[1]);
		block_size = atoi(argv[2]);
	}	
	A = getA15(N);
	ipiv = (int *)malloc(N * sizeof(int));
	int n = N;
	int info = -10;

	gpu_based_matrix_t gpu_matrix;

	checkCudaErrors(cudaGetDeviceCount(&(gpu_matrix.gpu_n)));
	gpu_matrix.order = N;
	gpu_matrix.block_size = block_size;

	cublasHandle_t* handle = (cublasHandle_t *)malloc(gpu_matrix.gpu_n * sizeof(cublasHandle_t));
	//cudaStream_t* stream = (cudaStream_t *)malloc(gpu_matrix.gpu_n * sizeof(cudaStream_t));
	
	if(initCublasHandles(handle, gpu_matrix.gpu_n))
	{
		printf("CUBLAS HANDLE INITIALIZATION FAILED\n");
		return EXIT_FAILURE;
	}

	/*if(initCudaStreams(stream, gpu_matrix.gpu_n))
	{
		printf("CUDA STREAM INITIALIZATION FAILED\n");
		return DEPT150_CUBLAS_ERROR;
	}*/
	
	if(copyMatrixToDevices(A, &gpu_matrix,handle ))
	{
		printf("INITIALIZATION OF DEVICE MATRIX FAILED");
		return EXIT_FAILURE;
	}
	/*printf("Before\n");
	print_matrix(N, N, A, N);
	printf("================\n");*/
	clock_t start_time, end_time;
	start_time = clock();
	double *b;
	int res = gaus_factorize_gpu_based_matrix(gpu_matrix, handle, /*stream,*/ ipiv);	
	end_time = clock();

	double *result_matrix = (double*)malloc(N*N*sizeof(double));
	copyMatrixFromDevices(result_matrix, gpu_matrix);

	printf("Start Checking result\n");
	double *proudct_LAndU = (double*)malloc(N*N*sizeof(double));
	bool hasErrors = false;
	for(int i = 0; i < N; i++)
	{
		for(int j =0; j < N; j++)
		{
			double sum = 0.;
			for(int k = 0; k < min(i, j); k++)
			{
				sum+= result_matrix[IDF2(i, k, N)] * result_matrix[IDF2(k, j, N)]; 
			}
			
			if (j>=i)
			{
				sum+= result_matrix[IDF2(i, j, N)];
			}else
			{
				sum+= result_matrix[IDF2(i, j, N)] * result_matrix[IDF2(j, j, N)]; 
			}

			proudct_LAndU[IDF2(i, j, N)] = sum;
		}
	}

	for(int i = 0; i < N; i++)
	{
		for(int j =0; j < N; j++)
		{
			if (fabs( proudct_LAndU[IDF2(i, j, N)] - A[IDF2(ipiv[i], j,N)] ) > 1e-6)
			{
				hasErrors = true;
				printf("For element [%i, %i] value of product %f initial value %f  Aindex %i\n", i, j, proudct_LAndU[IDF2(i, j, N)], A[IDF2(ipiv[i], j,N)], IDF2(ipiv[i], j,N)); 
			}
		}
	}

	if(!hasErrors)
	{
		printf("Comparing result product success\n");
	}
	
	//printf("Matrix after factorize \n");
	//print_matrix(N, N, result_matrix, N);

	/*printf("L and U product \n");
	print_matrix(N, N, proudct_LAndU, N);*/

	/*printf("IPIV \n");
	for(int i = 0 ; i < N; i++)
	{
		printf("%i ", ipiv[i]);
	}*/


	printf("\nTime : %.5f\n",(float)(end_time - start_time)/CLOCKS_PER_SEC);
	
	free(A);
	free(result_matrix);

	cublasStatus_t stat;
	for(int i = 0; i < gpu_matrix.gpu_n; i++)
	{
		checkCudaErrors(cudaSetDevice(i));
		checkCudaErrors(cudaFree(gpu_matrix.gpu_matrix_blocks[i]));
	}

	destroyCublasHandles(handle, gpu_matrix.gpu_n);

	free(gpu_matrix.gpu_matrix_blocks);
}
